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Robust empirical constitutive laws for granular materials in air or in a viscous fluid have been expressed in 
terms of timescales based on the dynamics of a single particle. However, some behaviours such as viscosity 
bifurcation or shear localization, observed also in foams, emulsions, and block copolymer cubic phases, seem 
to involve other micro-timescales which may be related to the dynamics of local particle reorganizations. In the 
present work, we consider a Tl process as an example of a rearrangement. Using the Soft dynamics simulation 
method introduced in the first paper of this series, we describe theoretically and numerically the motion of four 
elastic spheres in a viscous fluid. Hydrodynamic interactions are described at the level of lubrication (Poiseuille 
squeezing and Couette shear flow) and the elastic deflection of the particle surface is modeled as Hertzian. The 
duration of the simulated T\ process can vary substantially as a consequence of minute changes in the initial 
separations, consistently with predictions. For the first time, a collective behaviour is thus found to depend on 
another parameter than the typical volume fraction in particles. 

PACS numbers: 02.70.Ns, 82.70.-y, 83.80.Iz 



I. INTRODUCTION 

Many materials are made of particles in a surrounding fluid. 
Among them foams, emulsions, granular matter, colloidal sus- 
pensions and micro gels are of daily use. A great deal of 
research revealed their complex behaviors including elastic, 
plastic and viscous characters (H [U, [H 0] . This complexity 
results from the wide range of particle properties and parti- 
cle interactions involved. Great hints to comprehensive The- 
ological models were obtained by considering t he dyn amics 
of a single particle. Thus emerged the time ^Jm/RP for a 
single grain of mass m, accelerated by the normal stress P 
(force PR*), to move over a distance comparable to its own 
size R JH |a 0], the time rj jP for a grain immersed in a fluid 
of viscosity r\ subjected to the same normal stress iH, and 
the relaxation time rj R 2 /a for a bubble or a droplet with sur- 
face tension a in a viscous fluid H Qi| EUl- The effective 
viscosity was expressed as an empirical function of these mi- 
croscopic timescales (fjj, thereby providing robust scaling ex- 
pressions for variou s proper ties of grains lfl2[fT3lfl4ll and bub- 
bles EHHBIilES]. 

Nevertheless, particulate materials exhibit some uncommon 
rheological properties which seem to involve other timescales. 
Oscillatory shear experiments [4], and more generally the 
delayed adaptation of the shear rate to a sudden change in 
the applied stress JH Q H reveal long 
internal relaxation processes. Other observations such as a 
critical shear rate below which no homogeneous flow ex- 
ists JH 0, d, 0, l28ll . or the coexistence of liquid and 
solid regions (shear localization, shear banding, cracks) in 
emulsions |2<1 [3(3], foams |[3ll [33l l34tl. wormlike mi- 
celles (H Um and granular materials ffjl Il4 also 
point to a complex internal dynamics. Usually, this internal 
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dynamics is qualitatively understood as the competition be- 
tween external solicitations that the particles ex peri ence and 
their ability to move within their neighborhood Il22ll . Such a 
mechanism is the core of the definition of the jamming tran- 
sition in glassy sy stems, which is a subject of intense de- 
bate JMISElElSH- Getting new insights into reorga- 
nization micro-timescale should therefore clarify the origin of 
such properties and should also provide useful hints to refine 
and generalize existing models of the material response. 

A common reorganization process is the separation of par- 
ticles while other particles approach and fill the void. When 
they involve four particles, these events, usually refered to as 
T\ processes when dealing with foams and emulsions, occur 
in deformed regions at a frequency proportional to the defor- 
mation rate (see for example |3ll^|4^[4^[4^[48[[49t[50i 
I5II HHl). They relax some stress and dissipate some energy. 
The relation between the duration of a T\ and the local stress is 
thus expected to affect the rheological behaviour of the mate- 
rial. For dry foams, the T\ dynamics has been shown to result 
from the surface tension and the surface viscosity |[53tl . The 
stretching ability of the interfaces avoids the need to squeeze 
violently the fluid between the approaching bubbles, and its 
viscous dissipation is thus negligible. By contrast, the T\ dy- 
namics is less well described in wet foams or other less con- 
centrated systems. A comprehensive description of their dy- 
namics requires a careful description of the particle interac- 
tion. For instance, visco-elastic and even adhesive properties 
of particles were shown to be important lf54l l55i I5r3l . 

In this paper, we show that squeezing the liquid between 
close particles (here in three dimensions) can give rise to long 
relaxation times. To this aim, we do not focus on a specific 
material which would include the interaction between solid 
grains, bubbles, droplets or colloidal particles. Rather, we 
address the ubiquitous situation of elastic-like particles in a 
Newtonian fluid. We consider a simple system of in-plane 
spheres undergoing a T\ process, as depicted on Fig. [TJ We 
discuss under which circumstances a T\ process should in- 
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deed occur and (if it does occur) the relative contribution to 
the dynamic of the normal approach and separation versus the 
tangential sliding of particles Ic33ll . In three dimensions, for a 
dry foam or a concentrated emulsion, such a T\ process with 
four topologically active particles can in fact be decomposed 
into two topologically simpler processes involving five parti- 
cles. However, whatever the exact process, the dynamics will 
still involve normal motions and tangential sliding (as well as 
rotation in general). Because normal motions are stronger, as 
we show below, we believe that no essential new phenomenon 
will emerge from other reorganization processes as compared 
to the time scale evidenced in the present work. 

This paper is the second of a series which presents the 
physics of materials made of close-packed elastic-like parti- 
cles immersed in a viscous fluid. In the first paper II57I1 . we 
focused on the normal separation of two particles, and we 
showed that the flow between their close surfaces interplays 
with the particle deformation in a non-trivial manner. As this 
feature was ignored so far in existing discrete element sim- 
ulations such as Molecular Dynamics f58ll (for elastic grains 
without a surrounding fluid) for Stokesian Dynamics 115911 (for 
non-deformable grains in a viscous fluid), we are introducing 
a new simulation method, named Soft-Dynamics, to account 
for it. In this paper, we include the tangential interaction, 
and we provide the main steps of the implementation of the 
Soft Dynamics method for the present context. This will serve 
as an introduction to the principle of larger scale simulations 
with this new method, which will include both particle ro- 
tation and boundary conditions, and which should constitute 
a promising tool for investigating the collective behaviors of 
many complex materials. 

As we shall see, the geometry addressed in the present pa- 
per, although rather symmetric (the centers of the three di- 
mensional particles are arranged within a plane, at the vertices 
of a losange), proves sufficiently rich to reveal how minute 
changes in the system configuration have an essential influ- 
ence upon its dynamics. 




FIG. 1 : Schematic representation of a T\ process with four in-plane 
spheres. Due to the applied forces, the group of four particles swap 
neighbours. Two particles separate while two other particles estab- 
lish contact. Meanwhile, the other particle pairs reorient, as shown 
by the evolution of angle 6 from about ^ to about ^ . 



II. MODELLING PARTICLE INTERACTIONS 

When addressing the question of a T\ process between elas- 
tic spheres in a viscous fluid, see Fig.Q] most of the interac- 
tions have already been described in the first paper lf57ll . The 
only new feature is particle sliding and, correspondingly, tan- 
gential forces. Hence, quantities such as viscous friction co- 
efficients or spring constants are now tensorial. We express 
these interactions in the present section. Let us recall that we 
deal with three dimensional particles. 



A. Pairwise interactions 

As discussed in detail in the first paper 03l> because we 
consider rather dense systems where each particle is close to 
several other particles (surface-to-surface gap much smaller 
than the particle size}, we simply discard long-range, many- 
body interactions I6411 . Furthermore, under such thin gap 
conditions, the interacting region between particles is much 
smaller than the particle size and the interactions between a 
particle and its neighbours are mostly independent from each 
other and can therefore be treated as a sum of pairwise inter- 
actions. 

Particle i is subjected to some force F ij by its neighbour- 
ing particles j, which can be decomposted into (i) the local 
pressure field in the fluid that results from a viscous lubrica- 
tion interaction and ( ii) a remote interaction (such as a damped 
electrostatic interaction, steric repulsion, van der Waals inter- 
actions, disjoining pressure, etc): 



B. Particle surface deflection 



(1) 



A fraction of the above force F ij exerted by particle j tran- 
sits through a small portion of the surface of particle i and 
deflects it elastically. 

In practice, the viscous component F Vls of the force is en- 
tirely transmitted by the surface of particle /'. The effect of the 
remote component F rem is more subtle. Electrostatic forces 
between surface charges act upon the surface and contribute 
entirely to the elastic deflection. By contrast, Van der Waals 
interactions also act directly within particle ;. However, most 
of such interactions occur within a depth comparable with 
the inter-particle gap, which is always much smaller than the 
depth of the region that is deformed elastically (see paragraph 
HIF]below). 

Hence, for simplicity, it is reasonable to assume that the 
total force F ij between both particles entirely contributes to 
the elastic surface deflection: 



pelf ~ p .. 



(2) 



The expression of F f" in terms of the corresponding surface 
deflection is discussed in paragraph lll Fl below. 
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C. Force balance for each particle 

The sum of all forces applied to particle i, both the external 
force F f 3 and the pairwise forces F y is equal to the mass m, 
times the acceleration. This is assumed to vanish due to the 
dominant effect of the fluid viscosity over inertia: 



Ff + £7?y=m i X i 







(3) 



where F f xt is an external force acting on grain i (such as grav- 
ity) and where the sum runs over the neighbours of particle i. 
In principle, there is another equation, similar to Eq. (01, for 
the torques applied to particle ;. But as mentioned earlier [63], 
this is not needed for the present symmetric T\ configuration 
such as that of Fig. Q] 

The Soft Dynamics method ll57tl simulates the dynamics of 
such a system, determined by the system of Eqs. (??) for all 
interactions and Eqs (O for all particles i. In the present work, 
for simplicity, we omit the remote interactions in Eq. (??) as 
we did before ll57ll . 

In order to specify the elastic and viscous forces, let us now 
describe the geometry and the kinematics of the interacting 
region between a pair of neighbouring particles. 



D. Contact geometry and kinematics 

Let i and j denote two interacting particles, as depicted on 
Fig. |2] As compared to the first paper, the positions of the 
particles centers are now vectors, labeled OXj and OXj, and 
X ij = OXj — OXj is the center-to-center vector. The deflec- 
tions of the particle surfaces are also vectors, labeled 8 \ and 

8 'j. Since all particles are identical and since the (lubrica- 
tion) forces are pairwise and act locally, facing deflections are 

symmetric: Sj + 8'j = 0. Thus, for simplicity, we shall use 

the total deflection 8 ij = Sj — 8'j for each pair of interact- 
ing particles. The unit vector normal to the contact can be 
expressed as 



X, 



(4) 



The gap h,j between both particle surfaces depends on both 
the center- to-center vector Xy and the total deflection 5 y: 



h i j = \X ij - 8ij\-2R, 



(5) 



Similarly, the relative velocity of the material points that con- 
stitute each particle surface, v s , involves the translation ve- 
locity of the particles (as already mentioned [63], the particles 
do not rotate in the present situation) and the evolution of the 
surface deflection: 



In order to specify viscous and elastic interactions, we will 
need to deal with projectors and tensors. We will use the 
symbol "■" for the tensor product (contraction of one coordi- 
nate index), and u will denote the transposed of vector u . 
Hence, ~u T ■ V" = ~v fT ■ ~u will be the scalar product of ~~u and 
~v ', and ~u ■ ~v fT their outer product, which is a tensor. In par- 
ticular, we will make use of tensor a defined as the projector 
onto the normal direction: 



a = n 



(7) 




FIG. 2: Model of interaction for two elastic spheres in a viscous 
fluid, (a) elastic deflection of the surfaces (in traction); (b) normal 
dissipation due to Poiseuille flow in the gap; (c) tangential dissipation 
due to the Couette flow. The force is transmitted from a particle to 
another through the fluid and through a possible remote force. Such 
a system behaves like a Maxwell fluid (a dashpot and a spring in 
series). The effective friction is a function of the gap h and of the 
size a of the surface through which the force is transmitted. 



E. Viscous force 

For a pair of close spheres, as discussed earlier J57ll . the 
fluid region that mediates most of the force between both par- 
ticles has a large aspect ratio, and the flow is essentially par- 
allel to the solid surfaces: the lubrication approximation can 
be used (see for example |60tl ). As before, the fluid inertia 
is negligible (low Reynolds numbers) and the viscous force 
F ' " acting on the surfaces depends linearly on their relative 
velocity v s - 



pvis 
'/. = 



37irja A 
Ttrja 2 



(8) 
(9) 

(10) 
(ID 



X, 



Su. 



(6) 



where the interparticle friction tensor Z has two components 
(normal and in-plane), expressed in terms of the unity tensor 
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1 and the projector a defined by Eq. (0. The normal viscous 
friction £ is related to the Poiseuille flow induced by squeez- 
ing or pulling lf57ll (see Fig. Wp), while the in-plane friction 
coefficient A reflects the tangential motion (sliding) between 
both particles, which generates a Couette (shear) flow in the 
gap (see Fig.[2fc). 



F. Elastic force 

Let us assume that the size a (discussed in the next para- 
graph) of the interacting region between particles i and j is 
known. Then, as before 115711 . the force depends linearly on 
the surface deflection, but this time the relation is tensorial: 



previously [57], this generates a Maxwellian contact dynam- 
ics through the combination of the elastic surface deflection 
and the viscous response of the fluid in the gap: it is possi- 
ble to move the center-to-center distance X ,- ; while keeping 

constant the deflection 8 ij, and vice-versa. But as compared 
to a classical Maxwell behaviour, the elastic element does al- 
ways behave linearly (Hertzian contact in the strong deflec- 
tion regime), and the viscous element does not have a constant 
value, as it depend on the geometry of the gap, see Eqs. (ISlfTTTi. 

The Soft-Dynamics method consists in calculating the rate 
of change of all center positions OXi and all gap deflections 

8 ij as a function of their current values, and integrating them 
over a small time step. 



F da = aE - 8, 

E = E(c n a + c'(l-a)) 



(12) 
(13) 



where the tensorial proportionality constant E is essentially 
the (scalar) Young modulus E, but incorporates geometrical 
constants on the order of unity c" and c' for the normal and 
tangential responses, respectively. 

The elastic response of bubbles and droplets were found to 
deviate from such a Hertz elasticity j54l p5tl , Although they 
have no bulk elasticity, the surface tension (7 confers them 
some elastic-like properties, and the elastic-like force mainly 
depends on the deflection 8, the size a of the interacting region 
and an effective Young modulus which scales like a/R. 



G. Size of the interacting region 



The size of the interacting region, again 15711 . depends ei- 
ther on the gap thickness h (Poiseuille regime) when the parti- 
cle surface is weakly deflected, or on the normal force (Hertz 
regime) when the particle surface can be c onsid ered planar. In 
the first case, it can be expressed as a m \l2Rh. In the second 
case, it is essentially independent of the tangential force ll6lll 
and can thus be expressed in terms of the normal deflection: 



a w y/R\8"\ = yR \n ■ 8 |. As explained earlier 15J, for 
the purpose of the Soft Dynamics method, we interpolate be- 
tween both behaviours of a in a simple manner: 



a(h,S") = y/R(2h + \S»\ 



(14) 



The choice of this interpolation is not physically supported, 
but it does not affect assymptotic behavior in both limits. 



A. Equations of motion 

The system satisfies one equation per interaction, namely 
Eq. (??), and one equation per particle, namely Eq. (0. We 
shall now see how it is possible to derive equations of motion. 

For this, we need to express the unknowns velocities 8 ij and 

X ij in terms of the current state of the system. 

From Eqs. (0, (0 and (fT4l . it appears that the size a of the 
interacting region can be expressed as a function of Xjj and 
8 ij. It then follows from Eqs. (0 and (Y2\ that the elastic 
force F e,a can also be expressed as a function of Xy and 8 y: 



* ela 



* ela 



(Xij, 8ij) 



(15) 



As a result, its time-derivative F f" can be expressed as a sum 

two terms: one of them is linear in X ij while the other is linear 

in 8 ij. The (tensorial) coefficient of each of these two terms 
is a function of the current system configuration, i.e., of all 

particle and gap variables OXi and 8 ij. Now, it follows from 
Eqs. (0, (01 and (??) that 8 ij is an affine function of X y : 

A 



p ela W rem ] 

ij ~ ij J 



(16) 



— / — y 

where Zy, F f" and F y depend on the current system con- 
figuration. Hence, F f" can be expressed as an affine function 
of Xa: 



' ela 



Gjj ■ X jj 



bij 



(17) 



III. METHOD OF THE SOFT-DYNAMICS SIMULATION 

The Soft-Dynamics method aims at simulating the time 
evolution of a system of elastic particles and in a viscous 
fluid, such as depicted in previous sections. Like usual dis- 
crete simulation methods, the motion of each particle center 
results from the force balance, Eq. (0. The specificity is that 
the interaction evolution results from the decomposition of 
the center-to-center distance given by Eq. (0. As illustrated 



where the coefficients Gy and b ij depend only on the cur- 
rent system configuration. The detailed calculation of these 
coefficients is provided in Appendix lAl 

From this, the time derivative of Eq. (0 yields a system of 
equations for the particle center velocities. The equation that 
corresponds to particle ; reads: 



X{Gij-& j -OX i )}=Yyb ij -F*< 



(18) 
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where the sums run over all neighbours of particle ;. 

— ► — y 

Note that because Gh = —Gjj and b h = — b jj, and if we 
assume that the sum of all external forces vanishes, 

£Ff* = 0, (19) 

then the sum of Eqs. ( TT8l ) for all particles i vanishes. In other 
words, these vector equations are not independent: one of 
them must be replaced, for instance, by the condition that the 
average particle velocity is zero: 

£QX,=0 (20) 

i 

Let us consider the system of Eq. d20b (or a similar one) to- 
gether with Eqs. (fT8b . taken for all particles i except one. This 
system of equations can be inverted to obtain the particle cen- 
ter velocities OXj. The gap velocities F fj a are then calculated 
from Eq. (T7\ . 



pair of particles (AD), implies some sliding of the diagonal 
pairs (see Fig.[TJ. 

At the early stages of the process, when 9 w 7t/3, the exter- 
nal forces A^ v and N y can be expressed in terms of the normal 
forces in the horizontal (A^) and diagonal (A^) pairs of parti- 
cles, and in terms of the sliding force Sj in the diagonal pairs: 

N x = N h + 2 l -N d -2^-S d (21) 

N y = 2^-N d + 2 l -S d (22) 

In fact, as we shall now see, the tangential force S d is much 
smaller than the normal forces. To show this, let us first notice 
that the tangential velocity is related to the angle 9 defined 
on Fig. [T] v, ~ — RQ, In the Poiseuille regime, the particles 
surfaces are weakly deflected and the horizontal and diagonal 
gaps are related to angle 9 through R+ j h% = (2R + h d ) cos 9 . 
Hence, the gap variations obey jh/, « h d cos 9 — 2R 9 sin 9, 



B. Choice of a numerical step 

Gaining the center velocity OX( requires to solve the linear 
system ( fT8b . Standard and efficient procedures are available 
to inverse it. We used a second order Newtonian scheme for 
the numerical integration of particle position and as well as 
deflections. A typical time in the problem is the Stokes time 
T taken by a single particle submitted to a typical force F to 
move over a distance R in a fluid with viscosity rj, see Eq. (1291 
below. The numerical time step is set to 10~ 3 in units of z 
for all simulations. Other numerical schemes, such as Runge 
Kutta method, should make simulations faster. Furthermore, a 
study of the optimal required time step will be necessary when 
dealing with significantly more than only four particles. 



IV. T\ DYNAMICS 

Let us now use the Soft-Dynamics method to simulate a 
single T\ process. The system is depicted on Fig. [3] initially, 
particles B and C are aligned horizontally, with a small gap ho, 
while particles A and D are aligned vertically. The diagonal 
gaps (between A and B, etc) have thickness ho too. 

A horizontal force A^ is applied on particles B and C while 
a vertical force Ny is applied on A and D. Various evolutions 
are possible depending on these two forces, which may or may 
not give rise to a T\ process (see Pig. [3j. Basically, a T\ occurs 
only if the interaction between particles B and C is tensile. The 
criterion for the occurrence of a T\ process will be derived 
below, as well as a scaling for its dynamics. The duration of a 
T\ will then be measured from the simulation. 



A. Theoretical predictions 

The T\ process, which consists in a separation of the hor- 
izontal pair of particles (BC) and an approach of the vertical 
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FIG. 3: Schematic evolution of four particles subjected to external 
forces. Force N x is horizontal and acts on particles B and D. Force 
N y is vertical and acts upon particles A and C. Both and N y can 
be either compressive (> 0) or tensile (< 0). Regimes (I) and (II) 
correspond to compressive forces. In regime (I), the configuration 
remains mostly unaltered. By contrast, a topological rearrangement 
(T\ process) occurs when N y > v3N x , which corresponds to region 
(II). When N x or N y is tensile, the four beads do not remain together, 
as can be seen from the time evolutions sketched for regimes (III)- 
(VI). On the whole, a T\ process always occurs when Ny > V$N X 
(regimes II, III and IV, white region). It is followed by particle sep- 
aration when A^.v is tensile (regimes III and IV). By contrast, no T\ 
process occurs when Ny < v3A^ (regimes I, V and VI, light grey 
region). 
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-/j/, f=a /i c / cos 6 + 2 v t sin 9 



(23) 



Let us now transform each term of the above equation by ex- 
pressing it as a function of the corresponding normal or tan- 
gential force by using the appropriate friction coefficient as 
defined by Eq. Q: 



2 C 



' C 2 



5rf 73 
' A 2 



(24) 



The relative magnitude of friction coefficients £ and A can be 
derived from Eqs. ( TTOb and (fTTT i: 



(25) 



A 2 V/j, 

where the size a of the interaction region is given by Eq. ( fT4b . 
We thus have C/A ~ ^//i in the Poiseuille regime and £/A f=a 
fl<S' ! /7z 2 in the Hertz regime. Hence, except for very large gaps 
h comparable to the particle size R, the normal friction is much 
larger than the sliding friction: £ ^> A. It follows that 



S d 



1 



■■h(N d -N h ) 



(26) 



2^C 

can be neglected in Eqs. ( T2T1422] |. Hence, the interaction force 
within the horizontal pair BC depends only on the applied 
forces: 



N h ^N x 



(27) 



This implies that, as pictured on Fig. [3] the gap will open 
and the T\ will proceed whenever N/, is tensile, i.e., when 
Ny N x \/3 (white region of the diagram). By contrast, the 
particles will not swap neighbours when iVy < N x v3 (light 
grey region). 

When Nf, is indeed tensile, we now wish to determine how 
long it takes for the horizontal pair of particles to separate. 

The dynamics of such a normal motion was detailed in 
Ref. lf57ll . Let us define the reduced force 



and the Stokes time 



T = 



ER 2 
6n7]R 2 



(28) 



(29) 



With the force N/, acting within the horizontal pair BC, the 
initial configuration (gap /20) corresponds to the Poiseuille 
regime if ho > hup and to the Hertz regime if Iiq < hup, where 



ihp 



= R K" 2 / 3 



(30) 



The corresponding rate of change of the gap ( ,571 1 can be ex- 
pressed as: 



h = 
h = 



h 
x 

xR 2 



(Poiseuille, h > hup) 



(Hertz, h < h HP ). 



(31) 
(32) 



Integrating these equations yields the typical time A required 
to achieve the separation of the horizontal pair BC of particles 
from an initial gap ho to a much larger gap hf w R: 



X In 



4 



hf 

ho 

(Poiseuille, Iiq > hup) 



Ti 2 



R z 

h 2 
"hp 



Tin 



h f 



4 R 2 
TIC'- T »T 

hi 



hHP, 
(Hertz, ho < hup) 



(33) 



(34) 



Once the gap hi, of the horizontal pair BC becomes compa- 
rable to R, the diagonal pairs such as AB slide rather quickly 
(since their A <C Q, and soon the gap /;,, of the vertical pair 
AD becomes significantly smaller than R. The time it then 
takes to reach the same value ho is again comparable to A. 

Hence, the order of magnitude given by Eqs. d33T[34l l for 
the time A is typically the expected order of magnitude for 
the duration of the entire T\ process. We will now test this 
prediction by comparing it with the simulation results. 



B. Result from simulations 

We implement the Soft Dynamics method to simulate a T\ 
process such as that depicted on Fig. [3] varying the two con- 



-3 



< 



trol parameters we pointed out above: the initial gap 10 
ho/R < 0.8 and the reduced force 10~ 4 < K < 0.1. For sim- 
plicity, there is no horizontal force (A^ v = 0). The reduced 
force given by Eq. d28l ) is then equal to K = \N y \/ER 2 \/3 and 
the Stokes time is T = 67Tv / 3t?^ 2 /|A^,| . 

Figure [4] displays the variations of several quantities in 
the course of a T\ process with a given set of parameters 
{ho = 10~ 2 R, K — 3.10 -3 ). In order to avoid discontinuities 
in the simulation, the force A^ v is increased from zero to its 
nominal value within a time T, and remains constant there- 
after. From a macroscopic point of view, for instance through 
the variation of the angle 9, the system seems to be almost 
blocked (9 « 7i/J) for a significant amount of time (f < 100t). 
It then starts moving to reach its final configuration {9 « 7T/6) 
where it remains thereafter (t > 250t). During the "blocked" 
phase, the applied force N y is transmitted through the diag- 
onal interaction such as AB, thereby inducing a tensile force 
Nj, « — Ny/ V3 in the horizontal pair BC. Hence, despite the 
overall "blocked" appearance of the system, the horizontal 
gap /i/, between particles B and C slowly increases from its 
initial value Iiq- Correspondingly, the horizontal friction de- 
creases. 

The fast moving period starts as soon as this friction is low 
enough. Particles B and C then separate quickly while parti- 
cles A and D in the vertical pair approach each other, thereby 
giving rise to sliding friction on the diagonal interactions. As 
particles A and D approach, the corresponding gap h v de- 
creases and the friction increases. This approach then slows 
down. Thus, although the system keeps moving, it appears to 
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reach a new "blocked" configuration, with no more sliding or 
horizontal traction, but only a vertical compression. 

As the system subjected to a constant force keeps moving, 
we need to arbitrarily define the end of the T\ process. Among 
various possible choices, we shall here consider that the T\ 
process is completed when the vertical interaction transmits 
most of the applied force (N v = Q.99N y ). The resulting du- 
ration of the T\ process is plotted on Fig. (other criteria 
would yield similar results). The first observation is that, for 
the range of initial gaps and particle stiffnesses we consider, 
the duration of the T\ is distributed over a wide range of time 
scales, roughly between 3t and 10 3 T. Next, we observe that 
these results match our theoretical predictions reasonably: 

• if the horizontal pair BC is in the Poiseuille regime, the 
T\ duration A scales like Tln(/?//jo). It thus depends 
on particle radius, on the applied force and on the fluid 
viscosity through T, as can be seen from Eq. 1291 . and 
slightly on the initial gap through the logarithmic factor. 
The T\ duration is then just a few times larger than the 
Stokes time t; 

• if pair BC is in the Hertz regime, the T\ duration 

A scales essentially like xk? , which implies a 

much stronger dependence on ho, and longer durations 
since the particles are soft. In this case, A can be much 
longer than the Stokes time x. 

Note that in the latter case, the separating pair of particles 
leave the Hertz regime and enter the Poiseuille regime in the 
late stages of separation (h > hup)- However, because the 
evolution is much slower in the Hertz regime, see Eqs. (|3T1 - 
|32| |. these late Poiseuille stages contribute very weakly to the 
overall T\ duration A. 

In summary, the numerical result for the duration of a T\ 
process presented on Fig.[5]are compatible with Eqs. d33l434b . 
They demonstrate that the duration of a T\ is hardly larger than 
the Stokes time T given by Eq. ( |29b as long as the surface de- 
flection is small compared to the inter-particle gap (Poiseuille 
regime) and thus depend mainly on the applied force, on the 
fluid viscosity and on the particle size. Remarkably, in the op- 
posite regime where the deflection is larger than the gap (Hertz 
regime), the T\ duration depends strongly on the interparticle 
gap and can reach very large values, as illustrated by Fig. [5] 

V. CONCLUSION: BEYOND VOLUME FRACTION 

In this paper, we studied one of the simplest reorganiza- 
tion processes for immersed, closed-packed, elastically de- 
formable particles in a simple geometry. We showed that the 
time needed for this process results principally from the vis- 
cous flow of the fluid into or out of the gap between pairs of 
almost contacting particles: it is always mostly driven by the 
normal approach or separation, while the role of tangential 
sliding is negligible. 

We also showed that the time needed can be very long when 
particles are close or soft (more explicitely, when the gap is 
much thinner than the particle surface deflection). This is the 



central result of the present study and, as we show below, it 
pleads towards going beyond the sole usual volume fraction 
to describe the state of a particulate material. 



A. Volume fraction and interparticle gap 

Let us consider four particles in a compact configuration 
such as that on Fig.QJ (angle 9 f=a 7r/3). More precisely, let 
us consider two variants of this configuration, with two dif- 
ferent values of the interparticle gap Iiq, say ho ~ 10~ 2 R and 
h ~ 10- 3 R. Let us now apply weak forces (say K ~ 10 2 ). 
In both situations, because the force is weak and the gaps 
are small, the center-to-center distances are almost identical. 
Hence, both situations cannot be distinguished at first sight. 

Yet it can be seen from Fig. [5] that the duration of the 7\ 
process will then differ substantially. 

Similarly, with a large, disordered assembly of grains, it 
is anticipated that there will exist different situations where 
the volume fraction is almost identical but where a change in 
the typical value of the interparticle gaps causes a dramatic 
alteration of the delay A after application of the stress for the 
system to set into appreciable motion. 

This conjecture will be tested in a future work, using simu- 
lations with a large number of particles. 



B. Dilatancy and permeation 

In even larger samples of granular materials in a compact 
state, it is anticipated that the need for some additional fluid 
to enable reorganization processes (a phenomenon called di- 
latancy and illustrated on Fig.Q~|i will become the main source 
of delay: fluid from loose or particle-free regions needs to 
permeate through the granular material which behaves as a 
porous medium [27|. Describing such a phenomenon requires 
introducing the liquid pressure, and is not included in our sim- 
ulation so far. 



C. Towards other materials 

As such, the present work applies to soft, plain, elastic par- 
ticles (such as elastomer beads or latex particles) immersed in 
a very viscous fluid. We showed that a possible physical ori- 
gin for a delay in the system response is the viscous flow in 
the thin gap between neighbouring particle surfaces. 

In other materials, however, other ingredients may also in- 
fluence this delay or even become dominant. For instance, for 
objects enclosed in fluid interfaces (vesicles, onions, bubbles, 
droplets, etc), phenomena such as Marangoni effects, surface 
viscosity, the dynamics of surfactant adsorption and Gibb's 
elasticity should play a role. Finer phenomena should also be 
considered, such as the hydrodynamics involved either near 
a moving "contact" line between two such objects or within 
Plateau borders. By contrast, for solid grains, very different 
phenomena may come into play, including solid friction. 
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FIG. 4: Time evolution of various quantities in the course of a T\ process, with parameters K = 31CP 3 , ho = 10 2 i? and N x = 0. The gap h 
(a), the normal and tangential forces (c) and the ratio \S n \/h (d) are plotted for the vertical (solid black lines, AD), the horizontal (dotted black 
lines, BC) and the diagonal (solid grey lines, AB etc) pairs of particles. On graph (c), the dotted grey line represents the tangential force Sj of 
a diagonal pair such as AB (which is zero for the vertical pair AD and horizontal pair BC). Graph (b) shows the angle 6 such as defined on 
Fig.Q] On graph (d), a pair of particles for which \8"\/h > 1 is in the Hertz regime. If \8"\/h < 1, it is in the Poiseuille regime. 




FIG. 5: (Color online) Typical duration A of a T\ process as a function of the initial gap /?o and of the dimensionless applied force K". The 
data points were obtained through the Soft-Dynamics simulation presented here. Blue open circles correspond to 7] s where the horizontal pair 
has remained in the Poiseuille regime during the entire process. Full red squares correspond to T\s, such as that represented on Fig. [4] whose 
horizontal (separating) pair has been in the Hertz regime for part of the time. The surface is that defined by the theoretical model for both 
regimes (Eqs. l33lj34l wifh h f = 2.5R). 



For each of these phenomena, a simplified yet realistic pair- 



wise interaction law will need to be expressed and can then be 
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included in the simulation rather easily. 



D. Perspectives 

The present study suggests that further investigations using 
the Soft-Dynamics method with larger systems (including par- 
ticle rotation as well as boundary conditions) should provide 
interesting results, not only with the present system of plain, 
elastic beads in a viscous fluid, but also with different types of 
particle interactions. By testing ideas such as the influence of 
the typical interparticle gap (or other quantities if the interac- 
tions are different), it should also provide hints for analytical 
modelling beyond the role of the particle volume fraction. 
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APPENDIX A: DYNAMICS OF PARTICLES 

In this Appendix, we deduce the particle dynamics, given 
by Eqs. ( fT8l . from the physical model of interactions and the 
mechanical equilibria described in Sec. [Ill] We start from the 
time derivative of the particle force balance given by Eq. ([3]): 



p ela . 
ij 



(Al) 



The second term involves the time derivative E = E{c n ~c t )6t 
of the contact stiffness expressed by Eq. ( fT3l l. Using ( 1A13) , it 
can be expressed as: 



E Ji 
a — • o 

2 



= -b 2 
aE (c„ — c t ) , 



a ■ 8 



(A6) 
(A7) 



Note that this term vanishes for c„ = c,. 

The third term involves the time derivative of a, which we 
express through its partial derivatives: a(8",h) = a§ n 8" +ai,h 
(we used the notation ag,, = and a/, = 4f ). Replac- 
ing 8" and h by their respective expressions in terms of X jj, 
Eqs. (IA14I) and ( IA15I ) lead to: 



E- 8 -> 
a— — =G 3 -Xij- b 3 



G 3 = 
~b 3 = 



2 

E-~8 



E 8 



a 8" ( ~Vs ~n — 8 



-a h h 



(A8) 
(A9) 

(A10) 



Finally, substituting the three results of Eqs. ( [A3) . dA6) 
and ( |A8) into Eq. dA2b yields: 



G-X,j= b 



• ela 



(All) 



where G = Gi + G3 and b = b 1 + b 2 + b 3 are two ex- 

— > — y 

plicit functions of 8 and X ij Then, for each particle i, sum- 
ming on its interacting particles j and using the force balance 
(0 yields the system of equation ( 118) . 



Let us express the above as a sum of ( i) F f which is sup- 
posed to be known, ( ii) terms that are linear in the particle cen- 
ter velocities, and ( Hi) another term that is explicitly known 

from the current state of the system, i.e., from X ,- ; and 8 . To 
this aim, using Eq. (TT21 . let us express F f" in terms of the 
partial derivative of F eIa (a,E, 8 ): 



pela 
'j 



E 8 +-E- 8 



E 8 



(A2) 



2. Preliminary differentiations 



According to the definition of the normal vector, n = 
i £ , and to that of the associated projector, a = ~n ■ ~n T , 



we obtain their time derivatives: 
n = (1 — a) 



(A12) 
(A13) 



1. Combination 

We can now easily express each terms of Eq. dA2b as a func- 
tion of X ij. The first term is directly given by (O: 
a 



-E- 8 =Gi-X, 

^ a _ 
G X = -E 

bi = G 1 -~v s 



(A3) 
(A4) 
(A5) 



as two an explicit functions of 8 and X , 7 . Indeed, according 
to Eq. ( TTol l. v s can be expressed as a function of the elastic 

force: v s = Z~ l 



ela p rem 



The evolution of the normal deflection 8" = 8 ■ n can be 
expressed as a function of X jj by using Eqs. (IA12) and ©: 



8 n 



8 ■ n 

X tj - n 



— >r _^ 
8 ■ n 

V/ • 7? 



— >T _^ 

8 ■ n 



(A14) 
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Finally, from Eq. ©, we deduce the expression of the gap 
evolution h as a function of X y: 

h=(x?j-~8)-~n+~v s -~n. (A15) 
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